#########################################################################
#########################################################################
############################## SETUP STUFF ##############################
#########################################################################
#########################################################################

###############
### Setup R ###
###############

### Load packages
library(foreign)
library(sandwich)
library(randomForest)
library(lmtest)
library(MASS)
library(sciplot)
library(lme4)
library(DAAG)
library(ordinal)
library(sciplot)
library(foreach)
library(doParallel)
library(effects)
library(DataCombine)
library(plyr)
library(doMC)
library(plm)

### Clear space
rm(list = ls())

### Setup multi-core processing
library(parallel)

# Calculate the number of cores
no_cores <- detectCores() - 1

# Initiate cluster
cl <- makeCluster(no_cores)

##################
### Setup Data ###
##################

### Load data. Data from Abouharb, Moyer, and Schmidt (2013), Cingranelli, Richards, and Clay (2009), and Keith (2012). See paper for data source details and descriptions.
setwd("")
dat <- read.dta("new_evidence.dta")

### Download data for Melton and Ginsburg (2014) from http://comparativeconstitutionsproject.org/download-data/
### Load data
mg <- read.dta("MeltonGinsburg2014_TSCS.dta")
head(mg)

### Merge datasets
dat <- merge(dat,mg,by=c("YEAR","COW"),all=TRUE)

### Download data for Chilton and Versteeg (2015) from https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/NC06GQ
### Load data
cr <- read.table("ConRightsData.tab", sep="\t", header=TRUE)
colnames(cr)[1] <- "YEAR"
colnames(cr)[3] <- "COW"
head(cr)

### Merge datasets
dat <- merge(dat,cr,by=c("YEAR","COW"),all=TRUE)

### Inspect data
colnames(dat)

### Create binary GWF measure
dat$GWFDEM <- as.numeric(0)
dat$GWFDEM[dat$gwf_nonautocracy=="democracy"] <- as.numeric(1)

### Fix coding error from Abouharb, Moyer, and Schmidt (2013)
summary(dat$emil)
dat$emil[dat$emil=="2"] <- 1

### Create several complete case subsets
attach(dat)
datcr <- na.omit(data.frame(UNREG,LJI,lag_lji,postsd,lag_postsd,latentdoirtmean,lag_latentdoirtmean,latentdoirtsd, lag_latentdoirtsd,gdppc,gdpchng,democracy,mean,sd,newdemocscore,logpop,trdgdp,popdnsty, iccpr,interstateintensity,civilwarintensity,YEAR,COW,lag_ASSN,ASSN,lag_FORMOV,FORMOV, lag_DOMMOV,DOMMOV,lag_SPEECH,SPEECH,lag_NEW_RELFRE,NEW_RELFRE,lag_WORKER,WORKER,ELECSD, lag_ELECSD,NEW_EMPINX,lag_NEW_EMPINX,emil,royal,lgdppc,gdpgrowth,lpopdnsty, lag_meda,GWFDEM))
dim(datcr)
range(datcr$YEAR)
length(unique(datcr$COW))

datmech <- na.omit(data.frame(UNREG,LJI,lag_lji,postsd,lag_postsd,latentdoirtmean,lag_latentdoirtmean,latentdoirtsd, lag_latentdoirtsd,gdppc,gdpchng,democracy,mean,sd,newdemocscore,logpop,trdgdp,popdnsty, iccpr,interstateintensity,civilwarintensity,YEAR,COW,lag_ASSN,ASSN,lag_FORMOV,FORMOV, lag_DOMMOV,DOMMOV,lag_SPEECH,SPEECH,lag_NEW_RELFRE,NEW_RELFRE,lag_WORKER,WORKER,ELECSD, lag_ELECSD,NEW_EMPINX,lag_NEW_EMPINX,emil,royal,lgdppc,gdpgrowth,lpopdnsty, lag_meda,GWFDEM,lag_djji_int,libmov,assem_assoc,tradeun_strike,press_express,relig,RightsTotal, polpart,polpart_year,tradeun_strike_year,assem_assoc_year,relig_year,press_express_year,libmov_year))
dim(datmech)
range(datmech$YEAR)
length(unique(datmech$COW))

datmechdj <- na.omit(data.frame(UNREG,LJI,lag_lji,postsd,lag_postsd,latentdoirtmean,lag_latentdoirtmean,latentdoirtsd, lag_latentdoirtsd,gdppc,gdpchng,democracy,mean,sd,newdemocscore,logpop,trdgdp,popdnsty, iccpr,interstateintensity,civilwarintensity,YEAR,COW,lag_ASSN,ASSN,lag_FORMOV,FORMOV, lag_DOMMOV,DOMMOV,lag_SPEECH,SPEECH,lag_NEW_RELFRE,NEW_RELFRE,lag_WORKER,WORKER,ELECSD, lag_ELECSD,NEW_EMPINX,lag_NEW_EMPINX,emil,royal,lgdppc,gdpgrowth,lpopdnsty, lag_meda,GWFDEM,lag_djji_int,libmov,assem_assoc,tradeun_strike,press_express,relig,RightsTotal, polpart,polpart_year,tradeun_strike_year,assem_assoc_year,relig_year,press_express_year,libmov_year,defactojudicialindependence))
dim(datmechdj)
range(datmechdj$YEAR)
length(unique(datmechdj$COW))

datff <- na.omit(data.frame(UNREG,LJI,lag_lji,postsd,lag_postsd,latentdoirtmean,lag_latentdoirtmean,latentdoirtsd,lag_latentdoirtsd,gdppc,gdpchng,democracy,mean,sd,newdemocscore,logpop,trdgdp,popdnsty,iccpr,interstateintensity,civilwarintensity,YEAR,COW,lag_ASSN,ASSN,lag_FORMOV,FORMOV,lag_DOMMOV,DOMMOV,lag_SPEECH,SPEECH,lag_NEW_RELFRE,NEW_RELFRE,lag_WORKER,WORKER,ELECSD,lag_ELECSD,NEW_EMPINX,lag_NEW_EMPINX,emil,royal,lgdppc,gdpgrowth,lpopdnsty,mediascore,lag_meda,GWFDEM,fourfree))
dim(datff)
range(datff$YEAR)
length(unique(datff$COW))

#######################################################################
#######################################################################
############################## MAIN TEXT ##############################
#######################################################################
#######################################################################

### Create Plots ###
####################

par(mfrow=c(1,1))
plot(datmech$lag_lji,datmech$NEW_EMPINX,xlab="De Facto Judicial Independence (lagged)",ylab="Empowerment Rights Index",bty='n',col="#ADD8E6",pch=19)
abline(0,14,col="#666666",lty=2,lwd=4)
abline(lm(datmech$NEW_EMPINX~datmech$lag_lji),lwd=4)

emp_rights <- aggregate(NEW_EMPINX~YEAR, data=datmech, FUN=function(x) c(mean=mean(x,trim=.1),  se=se(x), count=length(x)))
jud_ind <- aggregate(lag_lji~YEAR, data=datmech, FUN=function(x) c(mean=mean(x,trim=.1),  se=se(x), count=length(x)))
old_jud_ind <- aggregate(defactojudicialindependence~YEAR, data=datmechdj, FUN=function(x) c(mean=mean(x,trim=.1),  se=se(x), count=length(x)))

cor(jud_ind$lag_lji,old_jud_ind$defactojudicialindependence)

emp_year <- emp_rights$YEAR
emp_upper <- emp_rights$NEW_EMPINX[,1]+1.64*emp_rights$NEW_EMPINX[,2]
emp_lower <- emp_rights$NEW_EMPINX[,1]-1.64*emp_rights$NEW_EMPINX[,2]
emp_rights_mean <- emp_rights$NEW_EMPINX[,1]
jud_year <- jud_ind$YEAR
jud_upper <- jud_ind$lag_lji[,1]+1.64*jud_ind$lag_lji[,2]
jud_lower <- jud_ind$lag_lji[,1]-1.64*jud_ind$lag_lji[,2]
jud_ind_mean <- jud_ind$lag_lji[,1]
old_jud_year <- old_jud_ind$YEAR
old_jud_upper <- old_jud_ind$defactojudicialindependence[,1]+1.64*old_jud_ind$defactojudicialindependence[,2]
old_jud_lower <- old_jud_ind$defactojudicialindependence[,1]-1.64*old_jud_ind$defactojudicialindependence[,2]
old_jud_ind_mean <- old_jud_ind$defactojudicialindependence[,1]

par(mfrow=c(1,2),
    oma = c(0,0,0,0),
    mar = c(5,5,2,2))
plot(emp_year,emp_rights_mean,xlab="Year",ylab="Empowerment Rights Index",bty='n',col="#FF8C77",type="l",lty=6,lwd="7",ylim=c(0,14))
lines(emp_year,emp_upper, lty=2, lwd=3)
lines(emp_year,emp_lower, lty=2, lwd=3)
plot(emp_year,jud_ind_mean,col="#ADD8E6",xlab="Year",ylab="De Facto Judicial Independence (Linzer and Staton 2015)",bty='n',type="l",lwd="7",ylim=c(0,1))
lines(emp_year, jud_upper, lty=2, lwd=3)
lines(emp_year, jud_lower, lty=2, lwd=3)

summary(datmech$lag_lji[datmech$YEAR==1982])
summary(datmech$lag_lji[datmech$YEAR==2008])
datmech[datmech$SPEECH==2 & datmech$YEAR==2008,]

plot(datmech$YEAR[datmech$COW==2],datmech$LJI[datmech$COW==2],xlim=c(1982,2008),ylim=c(0,1),xlab="Year",ylab="De Facto Judicial Independence",bty="n",type="l",lwd=6,col="#7fc97f") 
lines(datmech$YEAR[datmech$COW==710],datmech$LJI[datmech$COW==710],lwd=6,lty=2,col="#beaed4") 
lines(datmech$YEAR[datmech$COW==140],datmech$LJI[datmech$COW==140],lwd=6,lty=3,col="#fdc086") 
lines(datmech$YEAR[datmech$COW==290],datmech$LJI[datmech$COW==290],lwd=6,lty=4,col="#bf5b17") 
lines(datmech$YEAR[datmech$COW==560],datmech$LJI[datmech$COW==560],lwd=6,lty=5,col="#386cb0") 
lines(datmech$YEAR[datmech$COW==640],datmech$LJI[datmech$COW==640],lwd=6,lty=6,col="#f0027f") 
legend(1982,.1,c("U.S.","China","Brazil","Poland","S. Africa","Turkey"),lwd=c(6,6,6,6,6,6),lty=c(1,2,3,4,5,6),col=c("#7fc97f","#beaed4","#fdc086","#bf5b17","#386cb0","#f0027f"),cex=.5,horiz=TRUE, bty='n')

par(mar=c(5,4,4,5)+.1)
plot(datmech$YEAR[datmech$COW==2],datmech$LJI[datmech$COW==290],xlim=c(1982,2008),ylim=c(0,1),xlab="Year",ylab="De Facto Judicial Independence",bty="n",type="l",lwd=6,col="#bf5b17") 
par(new=TRUE)
plot(datmech$YEAR[datmech$COW==2],datmech$NEW_EMPINX[datmech$COW==290],axes=FALSE,ylim=c(0,14),xaxt="n",yaxt="n",xlab="",ylab="",bty="n",type="l",lty=3,lwd=6,col="#bf5b17")
axis(4)
mtext("Empowerment Rights Index",side=4,line=3)
legend(1991,4,c("De Facto Judicial Independence","Empowerment Rights Index"),lwd=c(6,6),lty=c(1,3),col=c("#bf5b17","#bf5b17"),cex=.7, bty='n')

par(mar=c(5,4,4,5)+.1)
plot(datmech$YEAR[datmech$COW==2],datmech$LJI[datmech$COW==640],xlim=c(1982,2008),ylim=c(0,1),xlab="Year",ylab="De Facto Judicial Independence",bty="n",type="l",lwd=6,col="#f0027f") 
par(new=TRUE)
plot(datmech$YEAR[datmech$COW==2],datmech$NEW_EMPINX[datmech$COW==640],axes=FALSE,ylim=c(0,14),xaxt="n",yaxt="n",xlab="",ylab="",bty="n",type="l",lty=3,lwd=6,col="#f0027f")
axis(4)
mtext("Empowerment Rights Index",side=4,line=3)
legend(1991,4,c("De Facto Judicial Independence","Empowerment Rights Index"),lwd=c(6,6),lty=c(1,3),col=c("#f0027f","#f0027f"),cex=.7, bty='n')

### Get Descriptives ###
########################

summary(datmech$NEW_EMPINX)
summary(datmech$lag_NEW_EMPINX)
summary(datmech$LJI)
summary(datmech$democracy)
summary(datmech$emil)
summary(datmech$royal)
summary(datmech$lgdppc)
summary(datmech$gdpgrowth)
summary(datmech$logpop)
summary(datmech$lpopdnsty)
summary(datmech$iccpr)
summary(datmech$interstateintensity)
summary(datmech$civilwarintensity)
summary(datmech$FORMOV)
summary(datmech$DOMMOV)
summary(datmech$SPEECH)
summary(datmech$ASSN)
summary(datmech$WORKER)
summary(datmech$NEW_RELFRE)
summary(datmech$ELECSD)
summary(datmech$libmov)
summary(datmech$assem_assoc)
summary(datmech$tradeun_strike)
summary(datmech$press_express)
summary(datmech$relig)
summary(datmech$polpart)

### Setup Data ###
##################

set.seed(12261991)
newdata <- list()
for(i in 1:1000){
  newdata[[i]] <- datcr
  newdata[[i]]$draw <- rnorm(cbind(rep(1,nrow(datcr))), mean=datcr$latentdoirtmean, sd=datcr$latentdoirtsd)
  newdata[[i]]$drawlag <- rnorm(cbind(rep(1,nrow(datcr))), mean=datcr$lag_latentdoirtmean, sd=datcr$lag_latentdoirtsd)
  newdata[[i]]$ljidraw <- rnorm(cbind(rep(1,nrow(datcr))), mean=datcr$lag_lji, sd=datcr$lag_postsd)
  newdata[[i]]$uds <- rnorm(cbind(rep(1,nrow(datcr))), mean=datcr$mean, sd=datcr$sd)
}

newdata2 <- list()
for(i in 1:1000){
  newdata2[[i]] <- datmech
  newdata2[[i]]$draw <- rnorm(cbind(rep(1,nrow(datmech))), mean=datmech$latentdoirtmean, sd=datmech$latentdoirtsd)
  newdata2[[i]]$drawlag <- rnorm(cbind(rep(1,nrow(datmech))), mean=datmech$lag_latentdoirtmean, sd=datmech$lag_latentdoirtsd)
  newdata2[[i]]$ljidraw <- rnorm(cbind(rep(1,nrow(datmech))), mean=datmech$lag_lji, sd=datmech$lag_postsd)
  newdata2[[i]]$uds <- rnorm(cbind(rep(1,nrow(datmech))), mean=datmech$mean, sd=datmech$sd)
  newdata2[[i]]$lji2 <- newdata[[i]]$ljidraw*newdata[[i]]$ljidraw
}

newdata3 <- list()
for(i in 1:1000){
  newdata3[[i]] <- datff
  newdata3[[i]]$draw <- rnorm(cbind(rep(1,nrow(datff))), mean=datff$latentdoirtmean, sd=datff$latentdoirtsd)
  newdata3[[i]]$drawlag <- rnorm(cbind(rep(1,nrow(datff))), mean=datff$lag_latentdoirtmean, sd=datff$lag_latentdoirtsd)
  newdata3[[i]]$ljidraw <- rnorm(cbind(rep(1,nrow(datff))), mean=datff$lag_lji, sd=datff$lag_postsd)
  newdata3[[i]]$uds <- rnorm(cbind(rep(1,nrow(datff))), mean=datff$mean, sd=datff$sd)
}

### Test with Composite CIRI Measure ###
########################################

xx <- c(1:14)

FML <- "as.factor(NEW_EMPINX) ~ lag_NEW_EMPINX + ljidraw + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(dat=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(15)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

### Plot model
var.names <- c("De Facto Judicial Independence (lagged)","Democracy","Military Regime","Monarchy","GDP Per Capita (logged)","GDP Growth (logged)",
               "Population (logged)", "Population Density (logged)","ICCPR Ratification","Interstate Conflict Intensity","Civil War Intensity")
var.names <- rev(var.names)
par(
  oma = c(0,0,0,0),
  mar = c(5,2,4,2)
)
plot(NULL,
     xlim = c(-2, 4),
     ylim = c(1, length(var.names)),
     axes = F, xlab = NA, ylab = NA)
abline(v = 0, lty = 3, lwd=2, col = "#aaaaaa")
est <- rev(model$beta[c(-1,-2)])
se <- rev(model$SE[c(-1,-2)])
for (i in 1:length(est)) {
  lines(c(est[i] + 1.96*se[i], est[i] - 1.96*se[i]), c(i, i), lwd = 15, col="#222222")
  points(est[i], i, pch = 19, cex = 1.1, col="white")
  points(est[i], i, pch = 19, cex = .7, col="black")
  text(est[i], i, var.names[i], xpd = T, cex = .8, pos = 3, offset=0.7)
}
axis(side = 1)
mtext(side = 1, "Estimated Coefficients", line = 2.5) 

lines(c(est[11] + 1.96*se[11], est[11] - 1.96*se[11]), c(11, 11), lwd = 15, col="#ADD8E6")
points(est[11], 11, pch = 19, cex = 1.1, col="white")
points(est[11], 11, pch = 19, cex = .7, col="black")

### Test with Individual CIRI Measures ###
##########################################

xx <- c(1:14)

FML <- "as.factor(FORMOV) ~ lag_FORMOV + ljidraw + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  or[,i] <- exp(tmp$coefficients)
  cil[,i] <- exp(confint(tmp)[,1])
  cih[,i] <- exp(confint(tmp)[,2])
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(15)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(DOMMOV) ~ lag_DOMMOV + ljidraw + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(15)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(SPEECH) ~ lag_SPEECH + ljidraw + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(15)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(ASSN) ~ lag_ASSN + ljidraw + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(15)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(WORKER) ~ lag_WORKER + ljidraw + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(15)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(NEW_RELFRE) ~ lag_NEW_RELFRE + ljidraw + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(15)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(ELECSD) ~ lag_ELECSD + ljidraw + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(15)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

### Plot models
var.names <- c("De Facto Judicial Independence (lagged) \n Model 2 (Foreign Movement)",
               "De Facto Judicial Independence (lagged) \n Model 3 (Domestic Movement)",
               "De Facto Judicial Independence (lagged) \n Model 4 (Speech)",
               "De Facto Judicial Independence (lagged) \n Model 5 (Assembly and Association)",
               "De Facto Judicial Independence (lagged) \n Model 6 (Workers' Rights)",
               "De Facto Judicial Independence (lagged) \n Model 7 (Religion)",
               "De Facto Judicial Independence (lagged) \n Model 8 (Electoral Self-Determination)")
betas <- c(2.270,2.250,3.790,3.738,2.634,3.109,4.003)
ses <- c(0.429,0.416,0.424,0.466,0.344,0.459,0.422)

var.names <- rev(var.names)

par(
  oma = c(0,0,0,0),
  mar = c(5,2,4,2)
)

plot(NULL,
     xlim = c(-1, 8),
     ylim = c(1, length(var.names)),
     axes = F, xlab = NA, ylab = NA)
abline(v = 0, lty = 3, lwd=2, col = "#aaaaaa")
est <- rev(betas)
se <- rev(ses)
for (i in 1:length(est)) {
  lines(c(est[i] + 1.96*se[i], est[i] - 1.96*se[i]), c(i, i), lwd = 15, col="#ADD8E6")
  points(est[i], i, pch = 19, cex = 1.1, col="white")
  points(est[i], i, pch = 19, cex = .7, col="black")
  text(est[i], i, var.names[i], xpd = T, cex = .8, pos = 3, offset=0.7)
}
axis(side = 1)
mtext(side = 1, "Estimated Coefficients", line = 2.5) 


### Test for Visibility Mechanism ###
#####################################

xx <- c(1:16)

FML <- "as.factor(FORMOV) ~ lag_FORMOV + ljidraw*SPEECH + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")    	   
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(17)]
  print(i)
}  # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(DOMMOV) ~ lag_DOMMOV + ljidraw*SPEECH + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")    	   
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(17)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(ASSN) ~ lag_ASSN + ljidraw*SPEECH + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")    	   
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(17)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(WORKER) ~ lag_WORKER + ljidraw*SPEECH + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")    	   
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(17)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean, na.rm=TRUE)
se.within <- apply(ses, 1, mean, na.rm=TRUE)
se.between <- apply(lms, 1, var, na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(NEW_RELFRE) ~ lag_NEW_RELFRE + ljidraw*SPEECH + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")    	   
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(17)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(ELECSD) ~ lag_ELECSD + ljidraw*SPEECH + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")    	   
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(17)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

### Test for Protection Mechanism ###
#####################################

xx <- c(1:16)

FML <- "as.factor(FORMOV) ~ lag_FORMOV + ljidraw*libmov + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(17)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(DOMMOV) ~ lag_DOMMOV + ljidraw*libmov + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(17)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(ASSN) ~ lag_ASSN + ljidraw*assem_assoc + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(17)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(WORKER) ~ lag_WORKER + ljidraw*tradeun_strike + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(17)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(SPEECH) ~ lag_SPEECH + ljidraw*press_express + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(17)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(NEW_RELFRE) ~ lag_NEW_RELFRE + ljidraw*relig + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(17)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "as.factor(ELECSD) ~ lag_ELECSD + ljidraw*polpart + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(17)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean,na.rm=TRUE)
se.within <- apply(ses, 1, mean,na.rm=TRUE)
se.between <- apply(lms, 1, var,na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

########################################################################
########################################################################
############################## APPENDICES ##############################
########################################################################
########################################################################

### Appendix A ###
##################

xx <- c(1:13)

FML <- "draw ~ drawlag + ljidraw + uds + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
for(i in 1:length(newdata)){
  tmp <- lmer(as.formula(FML), data=newdata[[i]])         
  lms[,i] <- fixef(tmp)
  ses[,i] <- sqrt(diag(vcov(tmp, useScale = FALSE)))
  print(i)
}
par.est <- apply(lms, 1, mean)
se.within <- apply(ses, 1, mean)
se.between <- apply(lms, 1, var)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results<- cbind(model$terms,model$beta,model$SE,model$beta/model$SE,2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "draw ~ drawlag + ljidraw + newdemocscore + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
for(i in 1:length(newdata)){
  tmp <- lmer(as.formula(FML), data=newdata[[i]])        
  lms[,i] <- fixef(tmp)
  ses[,i] <- sqrt(diag(vcov(tmp, useScale = FALSE)))
  print(i)
}
par.est <- apply(lms, 1, mean)
se.within <- apply(ses, 1, mean)
se.between <- apply(lms, 1, var)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results<- cbind(model$terms,model$beta,model$SE,model$beta/model$SE,2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

FML <- "draw ~ drawlag + ljidraw + GWFDEM + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
for(i in 1:length(newdata)){
  tmp <- lmer(as.formula(FML), data=newdata[[i]])        
  lms[,i] <- fixef(tmp)
  ses[,i] <- sqrt(diag(vcov(tmp, useScale = FALSE)))
  print(i)
}
par.est <- apply(lms, 1, mean)
se.within <- apply(ses, 1, mean)
se.between <- apply(lms, 1, var)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results<- cbind(model$terms,model$beta,model$SE,model$beta/model$SE,2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

par(
  oma = c(0,0,0,0),
  mar = c(5,2,4,2)
)

plot(NULL,
     xlim = c(-.1, 1),
     ylim = c(1, 4),
     axes = F, xlab = NA, ylab = NA)
abline(v = 0, lty = 3, lwd=2, col = "#aaaaaa")

lines(c(0.317,0.587), c(4, 4), lwd = 15, col="#ADD8E6")
points(c(0.452), 4, pch = 19, cex = 1.1, col="white")
points(c(0.452), 4, pch = 19, cex = .7, col="black")
text(c(0.452), 4, "De Facto Judicial Independence \n UDS Measure", xpd = T, cex = .8, pos = 3, offset=0.7)

lines(c(0.430,0.684), c(3, 3), lwd = 15, col="#ADD8E6")
points(c(0.557), 3, pch = 19, cex = 1.1, col="white")
points(c(0.557), 3, pch = 19, cex = .7, col="black")
text(c(0.557), 3, "De Facto Judicial Independence \n Demoracy-Dictatorship Measure", xpd = T, cex = .8, pos = 3, offset=0.7)

lines(c(0.395,0.653), c(2, 2), lwd = 15, col="#ADD8E6")
points(c(0.524), 2, pch = 19, cex = 1.1, col="white")
points(c(0.524), 2, pch = 19, cex = .7, col="black")
text(c(0.524), 2, "De Facto Judicial Independence \n Autocratic Regimes Measure", xpd = T, cex = .8, pos = 3, offset=0.7)

lines(c(0.467,0.721), c(1, 1), lwd = 15, col="#ADD8E6")
points(c(0.594), 1, pch = 19, cex = 1.1, col="white")
points(c(0.594), 1, pch = 19, cex = .7, col="black")
text(c(0.594), 1, "De Facto Judicial Independence \n Polity Measure", xpd = T, cex = .8, pos = 3, offset=0.7)

axis(side = 1)
mtext(side = 1, "Estimated Coefficients", line = 2.5) 

### Appendix B ###
##################

x <- c(cor(datmechdj$LJI[datmechdj$YEAR=="1982"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1982"]),cor(datmechdj$LJI[datmechdj$YEAR=="1983"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1983"]),cor(datmechdj$LJI[datmechdj$YEAR=="1984"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1984"]),cor(datmechdj$LJI[datmechdj$YEAR=="1985"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1985"]),cor(datmechdj$LJI[datmechdj$YEAR=="1986"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1986"]),cor(datmechdj$LJI[datmechdj$YEAR=="1987"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1987"]),cor(datmechdj$LJI[datmechdj$YEAR=="1988"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1988"]),cor(datmechdj$LJI[datmechdj$YEAR=="1989"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1989"]),cor(datmechdj$LJI[datmechdj$YEAR=="1990"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1990"]),cor(datmechdj$LJI[datmechdj$YEAR=="1991"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1991"]),cor(datmechdj$LJI[datmechdj$YEAR=="1992"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1992"]),cor(datmechdj$LJI[datmechdj$YEAR=="1993"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1993"]),cor(datmechdj$LJI[datmechdj$YEAR=="1994"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1994"]),cor(datmechdj$LJI[datmechdj$YEAR=="1995"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1995"]),cor(datmechdj$LJI[datmechdj$YEAR=="1996"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1996"]),cor(datmechdj$LJI[datmechdj$YEAR=="1997"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1997"]),cor(datmechdj$LJI[datmechdj$YEAR=="1998"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1998"]),cor(datmechdj$LJI[datmechdj$YEAR=="1999"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="1999"]),cor(datmechdj$LJI[datmechdj$YEAR=="2000"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="2000"]),cor(datmechdj$LJI[datmechdj$YEAR=="2001"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="2001"]),cor(datmechdj$LJI[datmechdj$YEAR=="2002"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="2002"]),cor(datmechdj$LJI[datmechdj$YEAR=="2003"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="2003"]),cor(datmechdj$LJI[datmechdj$YEAR=="2004"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="2004"]),cor(datmechdj$LJI[datmechdj$YEAR=="2005"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="2005"]),cor(datmechdj$LJI[datmechdj$YEAR=="2006"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="2006"]),cor(datmechdj$LJI[datmechdj$YEAR=="2007"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="2007"]),cor(datmechdj$LJI[datmechdj$YEAR=="2008"],datmechdj$defactojudicialindependence[datmechdj$YEAR=="2008"]))

par(mfrow=c(1,1),
    oma = c(0,0,0,0),
    mar = c(6,5,2,2))
plot(emp_year,x,bty='n',pch=19,col="#7A9F35",type="l",lwd="7",ylim=c(min(.65),max(.90)),xlab="Year",ylab="Correlation between the Keith (2012) and \n Linzer and Staton (2015) De Facto Judicial Independence Measures")

### Appendix C ###
##################

par(mfrow=c(1,2),
    oma = c(0,0,0,0),
    mar = c(5,5,2,2))
plot(emp_year,old_jud_ind_mean,col="#999999",xlab="Year",ylab="De Facto Judicial Independence (Keith 2012)",bty='n',type="l",lwd="7",ylim=c(0,2))
lines(emp_year, old_jud_upper, lty=2, lwd=3)
lines(emp_year, old_jud_lower, lty=2, lwd=3)
plot(emp_year,jud_ind_mean,col="#ADD8E6",xlab="Year",ylab="De Facto Judicial Independence (Linzer and Staton 2015)",bty='n',type="l",lwd="7",ylim=c(0,1))
lines(emp_year, jud_upper, lty=2, lwd=3)
lines(emp_year, jud_lower, lty=2, lwd=3)

### Appendix D ###
##################

xx <- c(1:13)

FML <- "draw ~ drawlag + ljidraw + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
for(i in 1:length(newdata)){
  tmp <- lmer(as.formula(FML), data=newdata[[i]])         
  lms[,i] <- fixef(tmp)
  ses[,i] <- sqrt(diag(vcov(tmp, useScale = FALSE)))
  print(i)
}
par.est <- apply(lms, 1, mean)
se.within <- apply(ses, 1, mean)
se.between <- apply(lms, 1, var)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results<- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

### Appendix E ###
##################

xx <- c(1:14)

FML <- "draw ~ drawlag + ljidraw + fourfree + uds + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata3))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata3))
for(i in 1:length(newdata3)){
  tmp <- lmer(as.formula(FML), data=newdata3[[i]])         
  lms[,i] <- fixef(tmp)
  ses[,i] <- sqrt(diag(vcov(tmp, useScale = FALSE)))
  print(i)
}
par.est <- apply(lms, 1, mean)
se.within <- apply(ses, 1, mean)
se.between <- apply(lms, 1, var)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata3))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results<- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

### Appendix F ###
##################

### Declare formulas and estimate OLS models
FORMULA <- "draw ~ drawlag"
FORMULA2 <- "draw ~ drawlag + ljidraw"
FORMULA3 <- "draw ~ drawlag + interstateintensity + civilwarintensity"
FORMULA4 <- "draw ~ drawlag + uds + emil + royal"
FORMULA5 <- "draw ~ drawlag + lgdppc + gdpgrowth + logpop + lpopdnsty"
FORMULA6 <- "draw ~ drawlag + uds + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity"
FORMULA7 <- "draw ~ drawlag + ljidraw + uds + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity"

### Declare functions
loss.rmse <- function(fit, df) cbind(as.numeric(predict(fit, newdata = df)), df$draw)
model.add <- function(df) lm(FORMULA, data=df)
model.add2 <- function(df) lm(FORMULA2, data=df)
model.add3 <- function(df) lm(FORMULA3, data=df)
model.add4 <- function(df) lm(FORMULA4, data=df)
model.add5 <- function(df) lm(FORMULA5, data=df)
model.add6 <- function(df) lm(FORMULA6, data=df)
model.add7 <- function(df) lm(FORMULA7, data=df)
validate.cv <- function(df, folds, resamples, model, loss) {
  lapply(1:resamples, function(x) {
    df$folds <- sample(rep(1:folds, length.out = nrow(df)))
    lapply(1:folds, function(test) {
      fit <- model(df[df$folds != test, ])
      loss(fit, df[df$folds == test, ])
    })})}

### Set values
SIMS <- 1000
cv1 <- NA
cv.add <- NA
cv2 <- NA
cv.add2 <- NA
cv3 <- NA
cv.add3 <- NA
cv4 <- NA
cv.add4 <- NA
cv5 <- NA
cv.add5 <- NA
cv6 <- NA
cv.add6 <- NA
cv7 <- NA
cv.add7 <- NA
newdata4 <- list()

for(i in 1:SIMS){
  
  newdata4[[i]] <- datcr
  newdata4[[i]]$draw <- rnorm(cbind(rep(1,nrow(datcr))), mean=datcr$latentdoirtmean, sd=datcr$latentdoirtsd)
  newdata4[[i]]$drawlag <- rnorm(cbind(rep(1,nrow(datcr))), mean=datcr$lag_latentdoirtmean, sd=datcr$lag_latentdoirtsd)
  newdata4[[i]]$ljidraw <- rnorm(cbind(rep(1,nrow(datcr))), mean=datcr$lag_lji, sd=datcr$lag_postsd)
  newdata4[[i]]$uds <- rnorm(cbind(rep(1,nrow(datcr))), mean=datcr$mean, sd=datcr$sd)
  
  cv.add <- validate.cv(newdata4[[i]], 10, 1, model.add, loss.rmse)
  cv.add2 <- validate.cv(newdata4[[i]], 10, 1, model.add2, loss.rmse)
  cv.add3 <- validate.cv(newdata4[[i]], 10, 1, model.add3, loss.rmse)
  cv.add4 <- validate.cv(newdata4[[i]], 10, 1, model.add4, loss.rmse)  
  cv.add5 <- validate.cv(newdata4[[i]], 10, 1, model.add5, loss.rmse)  
  cv.add6 <- validate.cv(newdata4[[i]], 10, 1, model.add6, loss.rmse)  
  cv.add7 <- validate.cv(newdata4[[i]], 10, 1, model.add7, loss.rmse)  
  
  value <- rbind(cv.add[[1]][[1]],cv.add[[1]][[2]],cv.add[[1]][[3]],cv.add[[1]][[4]],cv.add[[1]][[5]],cv.add[[1]][[6]],cv.add[[1]][[7]],cv.add[[1]][[8]],cv.add[[1]][[9]],cv.add[[1]][[10]])
  cv1[i] <- sqrt(mean((value[,1] - value[,2])^2))
  value2 <- rbind(cv.add2[[1]][[1]],cv.add2[[1]][[2]],cv.add2[[1]][[3]],cv.add2[[1]][[4]],cv.add2[[1]][[5]],cv.add2[[1]][[6]],cv.add2[[1]][[7]],cv.add2[[1]][[8]],cv.add2[[1]][[9]],cv.add2[[1]][[10]])
  cv2[i] <- sqrt(mean((value2[,1] - value2[,2])^2))
  value3 <- rbind(cv.add3[[1]][[1]],cv.add3[[1]][[2]],cv.add3[[1]][[3]],cv.add3[[1]][[4]],cv.add3[[1]][[5]],cv.add3[[1]][[6]],cv.add3[[1]][[7]],cv.add3[[1]][[8]],cv.add3[[1]][[9]],cv.add3[[1]][[10]])
  cv3[i] <- sqrt(mean((value3[,1] - value3[,2])^2))
  value4 <- rbind(cv.add4[[1]][[1]],cv.add4[[1]][[2]],cv.add4[[1]][[3]],cv.add4[[1]][[4]],cv.add4[[1]][[5]],cv.add4[[1]][[6]],cv.add4[[1]][[7]],cv.add4[[1]][[8]],cv.add4[[1]][[9]],cv.add4[[1]][[10]])
  cv4[i] <- sqrt(mean((value4[,1] - value4[,2])^2))
  value5 <- rbind(cv.add5[[1]][[1]],cv.add5[[1]][[2]],cv.add5[[1]][[3]],cv.add5[[1]][[4]],cv.add5[[1]][[5]],cv.add5[[1]][[6]],cv.add5[[1]][[7]],cv.add5[[1]][[8]],cv.add5[[1]][[9]],cv.add5[[1]][[10]])
  cv5[i] <- sqrt(mean((value5[,1] - value5[,2])^2))
  value6 <- rbind(cv.add6[[1]][[1]],cv.add6[[1]][[2]],cv.add6[[1]][[3]],cv.add6[[1]][[4]],cv.add6[[1]][[5]],cv.add6[[1]][[6]],cv.add6[[1]][[7]],cv.add6[[1]][[8]],cv.add6[[1]][[9]],cv.add6[[1]][[10]])
  cv6[i] <- sqrt(mean((value6[,1] - value6[,2])^2))
  value7 <- rbind(cv.add7[[1]][[1]],cv.add7[[1]][[2]],cv.add7[[1]][[3]],cv.add7[[1]][[4]],cv.add7[[1]][[5]],cv.add7[[1]][[6]],cv.add7[[1]][[7]],cv.add7[[1]][[8]],cv.add7[[1]][[9]],cv.add7[[1]][[10]])
  cv7[i] <- sqrt(mean((value7[,1] - value7[,2])^2))
}

y <- c(((cv2-cv1)/cv1),((cv3-cv1)/cv1),((cv4-cv1)/cv1),((cv5-cv1)/cv1),((cv6-cv1)/cv1),((cv7-cv1)/cv1))
a <- c(rep(c("CV Model 1 \n De Facto Judicial \n Independence \n Only"),1000),rep(c("CV Model 2 \n Conflict \n Only \n"),1000),rep(c("CV Model 3 \n Regime Types \n Only \n "),1000),
       rep(c("CV Model 4 \n Socioeconomic \n Only \n"),1000),rep(c("CV Model 5 \n Controls \n \n"),1000),rep(c("CV Model 6 \n Full Model \n \n"),1000))
d <- cbind(as.numeric(y),a)

### Plot reduction in mean square error across models
par(
  oma = c(0,0,0,0),
  mar = c(5,4,2,1)
)
bargraph.CI(d[,2],as.numeric(d[,1]), main="", col=c("#add8e6","#eeeeee","#eeeeee","#eeeeee","#eeeeee","#add8e6"), font.main=1,
            border=NA, family="Helvetica", err.lty=3, ylim=c(min(as.numeric(d[,1])),0),err.width=.05,xlab="Models",ylab="Proportional Reduction in Mean Square Error")

### Appendix G ###
##################

xx <- c(1:12)

set.seed(12261991)
FML <- "draw ~ drawlag + ljidraw + uds + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity"
mse <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
np <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
for(i in 1:length(newdata)){
  tmp <- randomForest(as.formula(FML), data=newdata[[i]], importance=TRUE, ntree=2000)     
  mse[,i] <- tmp$importance[,1]
  np[,i] <- tmp$importance[,2]
  print(i)
}
mse.combined <- apply(mse, 1, mean)
np.combined <- apply(np, 1, mean)
results <- cbind(xx, mse.combined, np.combined)
results

FML <- "as.factor(FORMOV) ~ lag_FORMOV + ljidraw + uds + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity"
mse <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
np <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
for(i in 1:length(newdata)){
  tmp <- randomForest(as.formula(FML), data=newdata[[i]], importance=TRUE, ntree=2000)     
  mse[,i] <- tmp$importance[,4]
  np[,i] <- tmp$importance[,5]
  print(i)
}
mse.combined <- apply(mse, 1, mean)
np.combined <- apply(np, 1, mean)
results <- cbind(xx, mse.combined, np.combined)
results

FML <- "as.factor(DOMMOV) ~ lag_DOMMOV + ljidraw + uds + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity"
mse <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
np <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
for(i in 1:length(newdata)){
  tmp <- randomForest(as.formula(FML), data=newdata[[i]], importance=TRUE, ntree=2000)     
  mse[,i] <- tmp$importance[,4]
  np[,i] <- tmp$importance[,5]
  print(i)
}
mse.combined <- apply(mse, 1, mean)
np.combined <- apply(np, 1, mean)
results <- cbind(xx, mse.combined, np.combined)
results

FML <- "as.factor(SPEECH) ~ lag_SPEECH + ljidraw + uds + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity"
mse <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
np <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
for(i in 1:length(newdata)){
  tmp <- randomForest(as.formula(FML), data=newdata[[i]], importance=TRUE, ntree=2000)     
  mse[,i] <- tmp$importance[,4]
  np[,i] <- tmp$importance[,5]
  print(i)
}
mse.combined <- apply(mse, 1, mean)
np.combined <- apply(np, 1, mean)
results <- cbind(xx, mse.combined, np.combined)
results

FML <- "as.factor(ASSN) ~ lag_ASSN + ljidraw + uds + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity"
mse <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
np <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
for(i in 1:length(newdata)){
  tmp <- randomForest(as.formula(FML), data=newdata[[i]], importance=TRUE, ntree=2000)     
  mse[,i] <- tmp$importance[,4]
  np[,i] <- tmp$importance[,5]
  print(i)
}
mse.combined <- apply(mse, 1, mean)
np.combined <- apply(np, 1, mean)
results <- cbind(xx, mse.combined, np.combined)
results

FML <- "as.factor(WORKER) ~ lag_WORKER + ljidraw + uds + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity"
mse <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
np <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
for(i in 1:length(newdata)){
  tmp <- randomForest(as.formula(FML), data=newdata[[i]], importance=TRUE, ntree=2000)     
  mse[,i] <- tmp$importance[,4]
  np[,i] <- tmp$importance[,5]
  print(i)
}
mse.combined <- apply(mse, 1, mean)
np.combined <- apply(np, 1, mean)
results <- cbind(xx, mse.combined, np.combined)
results

FML <- "as.factor(NEW_RELFRE) ~ lag_NEW_RELFRE + ljidraw + uds + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity"
mse <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
np <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
for(i in 1:length(newdata)){
  tmp <- randomForest(as.formula(FML), data=newdata[[i]], importance=TRUE, ntree=2000)     
  mse[,i] <- tmp$importance[,4]
  np[,i] <- tmp$importance[,5]
  print(i)
}
mse.combined <- apply(mse, 1, mean)
np.combined <- apply(np, 1, mean)
results <- cbind(xx, mse.combined, np.combined)
results

FML <- "as.factor(ELECSD) ~ lag_ELECSD + ljidraw + uds + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity"
mse <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
np <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
for(i in 1:length(newdata)){
  tmp <- randomForest(as.formula(FML), data=newdata[[i]], importance=TRUE, ntree=2000)     
  mse[,i] <- tmp$importance[,4]
  np[,i] <- tmp$importance[,5]
  print(i)
}
mse.combined <- apply(mse, 1, mean)
np.combined <- apply(np, 1, mean)
results <- cbind(xx, mse.combined, np.combined)
results

#########################################################################
#########################################################################
############################## OTHER TESTS ##############################
#########################################################################
#########################################################################

### Test for Quadratic Relationship ###
#######################################

xx <- c(1:16)

FML <- "as.factor(NEW_EMPINX) ~ lag_NEW_EMPINX + ljidraw + lji2 + democracy + emil + royal +  lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|COW)"
lms <- matrix(dat=NA, nrow=(length(xx)), ncol=length(newdata2))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata2))
for(i in 1:length(newdata2)){
  tmp <- clmm(as.formula(FML), data=newdata2[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(17)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean, na.rm=TRUE)
se.within <- apply(ses, 1, mean, na.rm=TRUE)
se.between <- apply(lms, 1, var, na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata2))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results <- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

### Test with Year RE ###
#########################

xx <- c(1:14)

FML <- "as.factor(NEW_EMPINX) ~ lag_NEW_EMPINX + ljidraw + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity + (1|YEAR)"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
for(i in 1:length(newdata)){
  tmp <- clmm(as.formula(FML), data=newdata[[i]], link = "logit", threshold = "equidistant")         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(15)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean, na.rm=TRUE)
se.within <- apply(ses, 1, mean, na.rm=TRUE)
se.between <- apply(lms, 1, var, na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results<- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

### Test with No RE ###
#######################

xx <- c(1:12)

FML <- "as.factor(NEW_EMPINX) ~ lag_NEW_EMPINX + ljidraw + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity"
lms <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
ses <- matrix(data=NA, nrow=(length(xx)), ncol=length(newdata))
for(i in 1:length(newdata)){
  tmp <- polr(as.formula(FML), data=newdata[[i]])         
  lms[,i] <- tmp$coefficients
  ses[,i] <- sqrt(diag(vcov(tmp)))[-c(13:26)]
  print(i)
} # If/when loop breaks, print 'i', skip 'i', and resume loop at 'i+1'.
par.est <- apply(lms, 1, mean, na.rm=TRUE)
se.within <- apply(ses, 1, mean, na.rm=TRUE)
se.between <- apply(lms, 1, var, na.rm=TRUE)
se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(newdata))))
model <- list("terms"=names(tmp$coefficients), "beta"=par.est, "SE"=se.est, "coefs"=lms)   
results<- cbind(model$terms, model$beta, model$SE, model$beta/model$SE, 2 * pnorm(abs(model$beta/model$SE), lower.tail = FALSE))
results

### Test with No RE and Using Clustered Standard Errors ###
###########################################################

### Not an elegant way of doing this with R that takes into account the uncertainty in the LJI point estimates

clx <- function(fm, dfcw, cluster) {
    library(sandwich)
    library(lmtest)
    library(zoo)
    M <- length(unique(cluster))
    N <- length(cluster)
    dfc <- (M/(M-1))*((N-1)/(N-fm$rank)) # anpassung der freiheitsgrade
    u <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum))
    vcovCL <-dfc * sandwich (fm, meat = crossprod(u)/N) * dfcw
    coeftest(fm, vcovCL)
}

FML <- "NEW_EMPINX ~ lag_NEW_EMPINX + lag_lji + democracy + emil + royal + lgdppc + gdpgrowth + logpop + lpopdnsty + iccpr + interstateintensity + civilwarintensity"
tmp <- lm(as.formula(FML), data=datcr)
clx(tmp, 1, datcr$COW)    